home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / isobuild / segment.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-02-11  |  15.3 KB  |  701 lines
  1.  
  2. /* segment.c                Brian Tierney, LBL
  3.  *
  4.  *  for use with the isobuild program
  5.  *
  6.  *  contains routines for segmentting 3D data, either by thresholding,
  7.  *   or by a 3D connectivity method
  8.  */
  9.  
  10. /* $Id: segment.c,v 1.6 1992/01/31 02:05:45 tierney Exp $ */
  11.  
  12. /* $Log: segment.c,v $
  13.  * Revision 1.6  1992/01/31  02:05:45  tierney
  14.  * *** empty log message ***
  15.  *
  16.  * Revision 1.5  1992/01/30  20:05:03  davidr
  17.  * prior to Brian's changes
  18.  *
  19.  * Revision 1.4  1992/01/10  01:59:31  davidr
  20.  * works with triserv now
  21.  *
  22.  * Revision 1.2  1991/12/19  01:42:20  davidr
  23.  * added RCS identification markers
  24.  * */
  25.  
  26. static char rcsid[] = "$Id: segment.c,v 1.6 1992/01/31 02:05:45 tierney Exp $" ;
  27.  
  28. #include "isobuild.h"
  29.  
  30. typedef struct s_item {        /* stack used to simulate recursion for
  31.                  * segmentation */
  32.     short     i, j, k;
  33. }         STACK_ITEM;
  34.  
  35. int       stack_size;
  36.  
  37. STACK_ITEM *stack = NULL;
  38. int       sp;            /* stack pointer */
  39.  
  40.  
  41. /* a 'grid' is used to perform the segmentation. The meaning of the
  42.  *  grid values are as follows:
  43.  *     1 : thresholded object
  44.  *     2 : edge of object  (if DO_BRIDGES is on )
  45.  *     3 : flood filled object
  46.  */
  47.  
  48. /* #define SEG_DEBUG */
  49.  
  50. /* useful globals */
  51. int       min_size;
  52.  
  53. /***************************************************************/
  54.  
  55. threshold_block(binfo, thresh)
  56.     BLOCK_INFO *binfo;
  57.     Data_type thresh;
  58. {
  59.     register int i, j;
  60.     int       x, y;
  61.  
  62.     y = binfo->yloc;
  63.     for (j = 0; j <= binfo->height; j++) {
  64.     x = binfo->xloc;
  65.     for (i = 0; i <= binfo->width; i++) {
  66.  
  67.         if (x < xdim && y < ydim) {
  68.         if (binfo->dslice[y][x] >= thresh)
  69.             binfo->grid[y][x] = 1;
  70.         }
  71.         x++;
  72.     }
  73.     y++;
  74.     }
  75. }
  76.  
  77. /***************************************************************/
  78. int
  79. seg_slice_thresh(thresh, data2d, grid2d, sx, sy, ex, ey, slice)
  80.     Data_type thresh;
  81.     Data_type **data2d;
  82.     Grid_type **grid2d;
  83.     int       sx, sy, ex, ey, slice;
  84. {
  85.     /* simple thresholding to segment data */
  86.  
  87.     register int i, j;
  88.     int       num_on = 0;
  89.  
  90. #ifdef SEG_DEBUG2
  91.     fprintf(stderr, "segmenting slice %d \n", slice);
  92. #endif
  93.  
  94.     for (j = sy; j < ey; j++)
  95.     for (i = sx; i < ex; i++) {
  96.         if (grid2d[j][i] != 255) {    /* not input mask */
  97.         if (data2d[j][i] >= thresh) {
  98.             grid2d[j][i] = 1;
  99.             num_on++;
  100.         } else {
  101.             grid2d[j][i] = 0;
  102.         }
  103.         }
  104.     }
  105.  
  106.     if (VERBOSE2)
  107.     fprintf(stderr, "%d voxels selected from slice %d \n", num_on, slice);
  108.  
  109. #ifdef SEG_DEBUG2
  110.     Status("after thresholding...");
  111.     show_grid_slice(grid2d, sx, sy, ex, ey);
  112. #endif
  113.  
  114.     return (num_on);
  115. }
  116.  
  117. /****************************************************************/
  118. int
  119. flood_fill_segmentation(thresh, bridge)
  120.     Data_type thresh;
  121.     int       bridge;
  122. {
  123.     int       rval, nvoxels, try_cnt,slice;
  124.     void      init_grid();
  125.     int       clear_background_grid();
  126.     void      clear_background_image(), locate_surfaces(), show_slices();
  127.  
  128.     init_grid(0);
  129.  
  130.     if (VERBOSE)
  131.     Status("Locating objects...");
  132.  
  133.     stack_size = identify_objects(thresh);    /* label objects '1' */
  134.  
  135.     if (SMOOTH_GRID) {
  136.     for (slice=SZ; slice<EZ; slice++)
  137.         smooth_slice(slice);
  138.     }
  139.  
  140.     fprintf(stderr, "%d voxels selected \n", stack_size);
  141.  
  142.     min_size = stack_size / 20;
  143.     if (min_size > 4000)
  144.     min_size = 4000;    /* dont let recursion go too deep */
  145.  
  146.     fprintf(stderr, "Looking for objects at least %d voxels in size \n",
  147.         min_size);
  148.  
  149.     alloc_stack(stack_size);    /* stack shouldn't be larger than the number
  150.                  * of pixels found by identify_objects */
  151.  
  152.     locate_surfaces();        /* label edges of object '2' */
  153.  
  154.     try_cnt = 0;
  155.     while ((rval = locate_object(bridge, try_cnt)) <= 0) {
  156.     if (rval < 0)
  157.         return (-1);
  158.     try_cnt++;
  159.     }
  160.  
  161.     if (VERBOSE)
  162.     Status("Zeroing non-objects...");
  163.  
  164.     nvoxels = clear_background_grid();    /* grid = 0 if grid != 3 */
  165.  
  166.     fprintf(stderr, "%d voxels selected after flood-fill segmentation \n",
  167.         nvoxels);
  168.  
  169.     cfree((char *) stack);
  170.  
  171. #ifdef SEG_DEBUG
  172.     show_slices();
  173. #endif
  174.  
  175.     return (thresh);
  176. }
  177.  
  178. /***************************************************************/
  179. int
  180. locate_object(bridge, try_cnt)
  181.     int       bridge, try_cnt;
  182. {
  183.     int       rval, pcnt;
  184.     static int sx = 1, sy = 1, sz = 1;
  185.     int       object_fill(), get_seed_guess1(), get_seed_guess2();
  186.     void      reset_grid();
  187.  
  188.     if (try_cnt == 0) {
  189.     rval = get_seed_guess1(&sx, &sy, &sz);
  190.     } else {
  191.     rval = get_seed_guess2(&sx, &sy, &sz);
  192.     }
  193.     if (rval < 0) {        /* didn't find an object */
  194.     Status("Object not found.");
  195.     return (-1);        /* give up */
  196.     }
  197.     fprintf(stderr, "Trying seed location: %d, %d, %d \n", sx, sy, sz);
  198.  
  199.     pcnt = object_fill(sx, sy, sz, bridge);
  200.  
  201.     if (pcnt < min_size) {
  202.     fprintf(stderr, "Object at location %d,%d,%d too small (%d pixels) \n",
  203.         sx, sy, sz, pcnt);
  204.     reset_grid(sx, sy, sz);    /* reset before doing flood fill */
  205.     if (sx < xmax)
  206.         sx++;
  207.     if (sy < ymax)
  208.         sy++;
  209.     if (sz < zmax)
  210.         sz++;
  211.     return (0);        /* try again */
  212.     } else {
  213.     fprintf(stderr, "Slice %d: Object at location: (%d,%d);  %d pixels \n",
  214.         sx, sy, sz, pcnt);
  215.     return (1);        /* success */
  216.     }
  217. }
  218.  
  219. /***************************************************************/
  220. void
  221. reset_grid(c, r, f)
  222.     int       c, r, f;
  223. {
  224. /* restore grid to original state before check_object_size routine
  225.  *  by setting all locations marked 3 back to 1
  226.  */
  227.  
  228.     sp = 0;            /* initialize stack pointer */
  229.  
  230.     push(-1, -1, -1);        /* null stack */
  231.     do {
  232.  
  233. start:
  234.     grid[f][r][c] = 1;
  235.  
  236.     if ((f < zmax) && (grid[f + 1][r][c] == 3)) {
  237.         push(c, r, f);
  238.         f++;
  239.         goto start;
  240.     }
  241.     if ((f > 0) && (grid[f - 1][r][c] == 3)) {
  242.         push(c, r, f);
  243.         f--;
  244.         goto start;
  245.     }
  246.     if ((r < ymax) && (grid[f][r + 1][c] == 3)) {
  247.         push(c, r, f);
  248.         r++;
  249.         goto start;
  250.     }
  251.     if ((r > 0) && (grid[f][r - 1][c] == 3)) {
  252.         push(c, r, f);
  253.         r--;
  254.         goto start;
  255.     }
  256.     if ((c < xmax) && (grid[f][r][c + 1] == 3)) {
  257.         push(c, r, f);
  258.         c++;
  259.         goto start;
  260.     }
  261.     if ((c > 0) && (grid[f][r][c - 1] == 3)) {
  262.         push(c, r, f);
  263.         c--;
  264.         goto start;
  265.     }
  266.     pop(&c, &r, &f);
  267.  
  268.     } while (f >= 0);        /* neg i indicates empty stack */
  269.  
  270.     if (sp != 0)
  271.     Error("stack not empty.");
  272.  
  273.     return;
  274. }
  275.  
  276. /***************************************************************/
  277. void
  278. show_slices()
  279. {                /* for debugging  */
  280.     register int c, r, f;
  281.  
  282.     for (f = 100; f < 101; f++) {
  283.     fprintf(stderr, "\n frame #: %d \n", f);
  284.     for (r = 0; r < ydim; r++) {
  285.         for (c = 0; c < xdim; c++)
  286.         fprintf(stderr, "%3d", grid[f][r][c]);
  287.         fprintf(stderr, "\n");
  288.     }
  289.     }
  290. }
  291.  
  292. /***************************************************************/
  293. int
  294. identify_objects(thresh_value)    /* labels objects with a '1' */
  295.     Data_type thresh_value;
  296. {
  297.     int       num_on = 0, check_val;
  298.  
  299.     long      nvoxels;
  300.     register int i;
  301.     register Grid_type *gptr;
  302.     register Data_type *dptr;
  303.  
  304.     gptr = **grid;
  305.     dptr = **data;
  306.  
  307.     nvoxels = zdim * xdim * ydim;
  308.  
  309.     for (i = 0; i < nvoxels; i++) {
  310.     check_val = (int) dptr[i];
  311.  
  312.     if (check_val < thresh_value)
  313.         gptr[i] = 0;
  314.     else {
  315.         gptr[i] = 1;
  316.         num_on++;
  317.     }
  318.     }
  319.  
  320.     return (num_on);
  321. }
  322.  
  323. /***************************************************************/
  324. void
  325. locate_surfaces()
  326. {                /* mark all surfaces as 2 */
  327.     register int c, r, f;
  328.  
  329.     /* label the surface points as 2 -- only checking 2-d slice */
  330.  
  331.     for (f = 0; f < zdim; f++)
  332.     for (r = 0; r < ydim; r++)
  333.         for (c = 0; c < xdim; c++)
  334.         if (grid[f][r][c] == 1) {
  335.             if ((r > 0 && grid[f][r - 1][c] == 0) ||
  336.             (r < ymax && grid[f][r + 1][c] == 0) ||
  337.             (c > 0 && grid[f][r][c - 1] == 0) ||
  338.             (c < xmax && grid[f][r][c + 1] == 0)) {
  339.             grid[f][r][c] = 2;
  340.             }
  341.         }
  342. }
  343.  
  344. /***************************************************************/
  345. int
  346. clear_background_grid()
  347.  /*
  348.   * at this point, all points in the desired object should be 3, so set any
  349.   * other points to 0
  350.   */
  351. {
  352.     register int i;
  353.     long      nvoxels;
  354.     register Grid_type *gptr;
  355.     int       cnt = 0;
  356.  
  357.     gptr = **grid;
  358.     nvoxels = zdim * xdim * ydim;
  359.  
  360.     for (i = 0; i < nvoxels; i++) {
  361.     if (gptr[i] == 3) {
  362.         gptr[i] = 1;
  363.         cnt++;
  364.     } else
  365.         gptr[i] = 0;
  366.     }
  367.     return (cnt);
  368. }
  369.  
  370. /***************************************************************/
  371. int
  372. mark_edge_neighbors(c, r, f, x_flag)    /* if neighbors are an edge, mark as
  373.                      * being part of object */
  374.     register int c, r, f, x_flag;
  375. {
  376.     int       cnt = 0;
  377.  
  378.     if ((r > 0) && (grid[f][r - 1][c] == 2)) {
  379.     grid[f][r - 1][c] = 3;
  380.     cnt++;
  381.     }
  382.     if ((r < ymax) && (grid[f][r + 1][c] == 2)) {
  383.     grid[f][r + 1][c] = 3;
  384.     cnt++;
  385.     }
  386.     if ((c > 0) && (grid[f][r][c - 1] == 2)) {
  387.     grid[f][r][c - 1] = 3;
  388.     cnt++;
  389.     }
  390.     if ((c < xmax) && (grid[f][r][c + 1] == 2)) {
  391.     grid[f][r][c + 1] = 3;
  392.     cnt++;
  393.     }
  394.     if (x_flag) {        /* also check x-slice */
  395.     if ((f > 0) && (grid[f - 1][r][c] == 2)) {
  396.         grid[f - 1][r][c] = 3;
  397.         cnt++;
  398.     }
  399.     if ((f < zmax) && (grid[f + 1][r][c] == 2)) {
  400.         grid[f + 1][r][c] = 3;
  401.         cnt++;
  402.     }
  403.     }
  404.     return (cnt);
  405. }
  406.  
  407. /***************************************************************/
  408.  
  409. int
  410. get_seed_guess1(sx, sy, sz)
  411.     int      *sx, *sy, *sz;
  412. {
  413. /* looks for a good seed value starting near the center of the image */
  414.     register int f, r, c;
  415.  
  416.     int       bx, by, bz;
  417.  
  418.     bx = (xdim / 2) - 2;    /* start near center of data set */
  419.     by = (ydim / 2) - 2;
  420.     bz = (zdim / 2);
  421.  
  422.     if (VERBOSE2)
  423.     fprintf(stderr, "Guess1: Looking for seed loc starting at: %d,%d,%d \n",
  424.         bx, by, bz);
  425.  
  426.     for (f = bz; f < zmax; f++)
  427.     for (r = by; r < ymax; r++)
  428.         for (c = bx; c < xmax; c++)
  429.         if (grid[f][r][c] == 1) {
  430.             if ((grid[f][r][c - 1] == 1) &&
  431.             (grid[f][r][c + 1] == 1) &&
  432.             (grid[f][r - 1][c] == 1) &&
  433.             (grid[f][r + 1][c] == 1)
  434.             && (grid[f][r - 1][c - 1] == 1) &&
  435.             (grid[f][r + 1][c + 1] == 1) &&
  436.             (grid[f][r + 1][c - 1] == 1) &&
  437.             (grid[f][r - 1][c + 1] == 1)) {
  438.             *sx = c;
  439.             *sy = r;
  440.             *sz = f;
  441.             return (0);
  442.             }
  443.         }
  444.     return (-1);
  445. }
  446.  
  447. /***************************************************************/
  448.  
  449. int
  450. get_seed_guess2(sx, sy, sz)
  451.     int      *sx, *sy, *sz;
  452. {
  453.     /*
  454.      * looks for a good seed value starting at the upper left corner of the
  455.      * image
  456.      */
  457.  
  458.     register int f, r, c;
  459.     int       bx, by, bz;
  460.  
  461.     bx = *sx;
  462.     by = *sy;
  463.     bz = *sz;
  464.  
  465.     if (VERBOSE2)
  466.     fprintf(stderr, "Guess2: Looking for seed loc starting at: %d,%d,%d \n",
  467.         bx, by, bz);
  468.  
  469.     for (f = bz; f < zmax; f++)
  470.     for (r = by; r < ymax; r++)
  471.         for (c = bx; c < zmax; c++)
  472.         if (grid[f][r][c] == 1) {
  473.             if ((grid[f][r][c - 1] == 1) &&
  474.             (grid[f][r][c + 1] == 1) &&
  475.             (grid[f][r - 1][c] == 1) &&
  476.             (grid[f][r + 1][c] == 1)
  477.             && (grid[f][r - 1][c - 1] == 1) &&
  478.             (grid[f][r + 1][c + 1] == 1) &&
  479.             (grid[f][r + 1][c - 1] == 1) &&
  480.             (grid[f][r - 1][c + 1] == 1)
  481.             ) {
  482.             *sx = c;
  483.             *sy = r;
  484.             *sz = f;
  485.             return (0);
  486.             }
  487.         }
  488.     return (-1);
  489. }
  490.  
  491.  
  492. /***************************************************************/
  493.  
  494. void
  495. init_grid(val)
  496.     int       val;
  497.  /* sets the entire grid to 'val' */
  498. {
  499.     int       nvoxels;
  500.  
  501.     nvoxels = ydim * xdim * zdim;
  502.     if (VERBOSE)
  503.     Status("Initializing grid.");
  504.  
  505. #ifdef OLD
  506.     register int i;
  507.     gptr = **grid;
  508.     register Grid_type *gptr;
  509.  
  510.     for (i = 0; i < nvoxels; i++)
  511.     gptr[i] = val;
  512. #else
  513.     memset((char *) grid[0][0], val, nvoxels);
  514. #endif
  515.  
  516. }
  517.  
  518. /***************************************************************/
  519. int
  520. count_diag_neighbor_edges(c, r, f, x_flag)
  521.     register int c, r, f, x_flag;    /* max value returned is 12 */
  522. {
  523.     /*
  524.      * NOTE: this routine only counts neighbors that are marked '2'. There
  525.      * will also be some neighbors marked 0, but we don't want to count those
  526.      * because when labeling the edges we only look at the 2-D slice.
  527.      */
  528.  
  529.     int       neighbors = 0;
  530.  
  531.     if ((r > 0) && (c > 0))
  532.     if (grid[f][r - 1][c - 1] == 2)
  533.         neighbors++;
  534.     if ((r > 0) && (c < xmax))
  535.     if (grid[f][r - 1][c + 1] == 2)
  536.         neighbors++;
  537.     if ((r < ymax) && (c > 0))
  538.     if (grid[f][r + 1][c - 1] == 2)
  539.         neighbors++;
  540.     if ((r < ymax) && (c < xmax))
  541.     if (grid[f][r + 1][c + 1] == 2)
  542.         neighbors++;
  543.  
  544.     if (x_flag) {        /* count x-slices too */
  545.     /* previous slice */
  546.     if ((f > 0) && (r > 0))
  547.         if (grid[f - 1][r - 1][c] == 2)
  548.         neighbors++;
  549.     if ((f > 0) && (c > 0))
  550.         if (grid[f - 1][r][c - 1] == 2)
  551.         neighbors++;
  552.     if ((f > 0) && (r < ymax))
  553.         if (grid[f - 1][r + 1][c] == 2)
  554.         neighbors++;
  555.     if ((f > 0) && (c < xmax))
  556.         if (grid[f - 1][r][c + 1] == 2)
  557.         neighbors++;
  558.     /* next slice */
  559.     if ((f < zmax) && (r > 0))
  560.         if ((grid[f + 1][r - 1][c] == 2))
  561.         neighbors++;
  562.     if ((f < zmax) && (c > 0))
  563.         if ((grid[f + 1][r][c - 1] == 2))
  564.         neighbors++;
  565.     if ((f < zmax) && (r < ymax))
  566.         if ((grid[f + 1][r + 1][c] == 2))
  567.         neighbors++;
  568.     if ((f < zmax) && (c < xmax))
  569.         if ((grid[f + 1][r][c + 1] == 2))
  570.         neighbors++;
  571.     }
  572.     return (neighbors);
  573. }
  574.  
  575. /***************************************************************/
  576. int
  577. object_fill(c, r, f, bridges)    /* label all points within the item 1 */
  578.     int       c, r, f, bridges;
  579.  
  580. /* this routine uses a stack to do a 3D recursive flood fill starting
  581.    at location (c,r,f)  */
  582. {
  583.     void      sum_pixel();
  584.  
  585.     int       cnt = 0, skip = 0;/* count number of points in object */
  586.     sp = 0;            /* initialize stack pointer */
  587.  
  588.     push(-1, -1, -1);        /* null stack */
  589.     do {
  590.  
  591. start:
  592.     grid[f][r][c] = 3;
  593.     cnt++;
  594.     if (bridges > 0 && cnt > 1) {    /* check strength of connecting
  595.                      * bridges */
  596.  
  597.         if (bridges == 1) {    /* 2D bridges */
  598.         cnt += mark_edge_neighbors(c, r, f, 0);
  599.         if (count_diag_neighbor_edges(c, r, f, 0) >= 2) {
  600.             skip++;
  601.             goto w_stop;
  602.         }
  603.         }
  604.         if (bridges == 2) {    /* 3D weak bridges */
  605.         cnt += mark_edge_neighbors(c, r, f, 1);
  606.         if (count_diag_neighbor_edges(c, r, f, 1) >= 4) {
  607.             skip++;
  608.             goto w_stop;
  609.         }
  610.         }
  611.         if (bridges == 3) {    /* 3D strong bridges */
  612.         cnt += mark_edge_neighbors(c, r, f, 1);
  613.         if (count_diag_neighbor_edges(c, r, f, 1) >= 2) {
  614.             skip++;
  615.             goto w_stop;
  616.         }
  617.         }
  618.     }
  619. #ifdef NEED?
  620.     else {
  621.         cnt += mark_edge_neighbors(c, r, f, 0);    /* include edges */
  622.     }
  623. #endif
  624.  
  625.     if ((f < zmax) && (grid[f + 1][r][c] == 1)) {
  626.         push(c, r, f);
  627.         f++;
  628.         goto start;
  629.     }
  630.     if ((f > 0) && (grid[f - 1][r][c] == 1)) {
  631.         push(c, r, f);
  632.         f--;
  633.         goto start;
  634.     }
  635.     if ((r < ymax) && (grid[f][r + 1][c] == 1)) {
  636.         push(c, r, f);
  637.         r++;
  638.         goto start;
  639.     }
  640.     if ((r > 0) && (grid[f][r - 1][c] == 1)) {
  641.         push(c, r, f);
  642.         r--;
  643.         goto start;
  644.     }
  645.     if ((c < xmax) && (grid[f][r][c + 1] == 1)) {
  646.         push(c, r, f);
  647.         c++;
  648.         goto start;
  649.     }
  650.     if ((c > 0) && (grid[f][r][c - 1] == 1)) {
  651.         push(c, r, f);
  652.         c--;
  653.         goto start;
  654.     }
  655. w_stop:
  656.     pop(&c, &r, &f);
  657.  
  658.     } while (f >= 0);        /* neg i indicates empty stack */
  659.  
  660.     if (sp != 0)
  661.     Error("stack not empty.");
  662.  
  663.     if (bridges)
  664.     fprintf(stderr, " fill stopped %d times due to weak bridges \n", skip);
  665.  
  666.     return (cnt);
  667. }
  668.  
  669. /***************************************************************/
  670. push(i, j, k)            /* add location to stack */
  671.     int       i, j, k;
  672. {
  673.     sp++;
  674.     if (sp >= stack_size) {
  675.     fprintf(stderr, "Stack was allocated %d slots \n", stack_size);
  676.     Error("Recursive stack overflow!!");
  677.     }
  678.     stack[sp].i = i;
  679.     stack[sp].j = j;
  680.     stack[sp].k = k;
  681. }
  682.  
  683. /***************************************************************/
  684. pop(i, j, k)            /* remove item from stack */
  685.     int      *i, *j, *k;
  686. {
  687.     *i = stack[sp].i;
  688.     *j = stack[sp].j;
  689.     *k = stack[sp].k;
  690.     sp--;
  691. }
  692.  
  693. /***************************************************************/
  694. alloc_stack(st_size)        /* allocate stack for non-recursive
  695.                  * flood-fill alg */
  696.     int       st_size;
  697. {
  698.     if ((stack = Calloc(st_size, STACK_ITEM)) == NULL)
  699.     perror("calloc: stack");
  700. }
  701.